The electromagnetic response in a layered vertical transverse isotropic medium: A new look at an old problem

We determined that the electromagnetic vertical transverse isotropic response in a layered earth can be obtained by solving two equivalent scalar equations, which were for the vertical electric field and for the vertical magnetic field, involving only a scalar global reflection coefficient. Besides the complete derivation of the full electromagnetic response, we also developed the corresponding computer code called EMmod, which models the full electromagnetic fields including internal multiples in the frequency-wavenumber domain and obtains the frequency-space domain solutions through a Hankel transformation by computing the Hankel integral using a 61-point Gauss-Kronrod integration routine. The code is able to model the 3D electromagnetic field in a 1D earth for diffusive methods such as controlled source electromagnetics as well as for wave methods such as ground penetrating radar. The user has complete freedom to place the source and the receivers in any layer. The modeling is illustrated with three examples, which aim to present the different capabilities of EMmod, while assessing its correctness. INTRODUCTION Solutions for waves and fields in layered earth models have been published in many forms and over many years. Here, the development of electromagnetic waves and fields is only briefly reviewed. The first to give the procedure for the scalar global reflection coefficient of a layered medium is Airy (1833) who analyzes the formation of Newton’s rings in a three-layered medium. For isotropic layered media with a source in the upper half-space, early treatments are given by Wait (1951, 1953) and for a horizontal electric dipole over a layered anisotropic half-space by Wait (1966). These results were obtained in the diffusive limit and for isotropic magnetic permeability. Redheffer (1961) introduces a linear fractional transformation, later known as the Redheffer star product, and uses it to formally solve the layered medium problem using the scattering matrix for electromagnetic waves. Kong (1972) solves the layered anisotropic half-space problem for arbitrary dipole sources in the upper half-space and considers both electric and magnetic vertical transverse isotropic (VTI) layers. He exploits the fact that the vertical components of the electric and magnetic fields are independent from each other. He still uses propagation matrices to derive the solutions, by invoking continuity of the horizontal components of the electric and magnetic field vectors. The problem with the propagation matrices method is that in the computational scheme, exponentially growing functions also occur, which makes it unusable for diffusive field models. Redheffer’s method uses the scattering matrix, which contains only exponentially damped functions and can be used for both wave and diffusive models. Tsang et al. (1975) solve the VTI-layered-earth problem in the high-frequency limit for microwave passive remote sensing of the earth. A solution for buried sources can be found in Ali and Mahmoud (1979). An excellent review on layered media models is given by Ursin (1983). Kwon and Wang (1986) solve the VTI-layeredearth problem by splitting transverse electric (TE) and transverse magnetic (TM) modes. Xiong (1989) uses vector potentials to solve for the electromagnetic fields in a stratified anisotropic medium. Michalski (2005) derives the Green’s tensor in layered VTI media using the transmission-line analogue. The scattering matrix formulation has been used recently to solve the forward problem for arbitrary anisotropic layered media with sources and receivers located in any layer of the model by Loseth and Ursin (2007). Zhong et al. (2008) and Wang et al. (2009) publish algorithms for induction logging tools. A code for forward and inverse modeling of isotropic low frequency electromagnetic data with complete freedom of placement of the source and receiver is published by Key (2009). Peer-reviewed code related to this article can be found at http://software.seg.org/2015/0001. Manuscript received by the Editor 8 November 2013; revised manuscript received 21 July 2014; published online 4 December 2014. Delft University of Technology, Department of Geotechnology, Delft, The Netherlands. E-mail: [email protected]; [email protected]; [email protected]. © 2014 Society of Exploration Geophysicists. All rights reserved. F1 GEOPHYSICS, VOL. 80, NO. 1 (JANUARY-FEBRUARY 2015); P. F1–F18, 6 FIGS. 10.1190/GEO2013-0411.1 D ow nl oa de d 04 /2 0/ 15 to 1 31 .1 80 .1 31 .2 42 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / In this paper, we show that the electromagnetic VTI response can be obtained by solving two equivalent scalar equations involving only a scalar global reflection coefficient. Unlike previous work featuring similar derivations, a corresponding computer code written in C and Fortran to model electromagnetic fields (EMmod) is attached to this paper. The combination of a simple derivation and a corresponding open-source code that is thoroughly tested (see the “Example” section), opens many possibilities for research and for benchmarking other codes. The code is able to model data for controlled-source electromagnetic (CSEM) methods in the diffusive limit as well as for ground-penetrating radar (GPR) methods for waves in dissipative media. The code allows placement of the source and the receiver anywhere in a stack of VTI layers for any source-receiver component. It outputs the electromagnetic field in the frequency-space domain. Whenever convenient, the subscript notation is used and the summation convention applies to repeated lowercase subscripts. Latin subscripts take on the values f1; 2; 3g, to represent vector directions along the three coordinates in space fx1; x2; x3g. Greek subscripts take on the values f1; 2g to denote vector directions along horizontal coordinates only. The unit matrix is denoted δkr defined as δkr 1⁄4 1 if k 1⁄4 r and δkr 1⁄4 0 otherwise. To denote vector products the Levi-Civita tensor of rank three is used, εklm 1⁄4 −1 when the sequence fklmg is an odd permutation of f123g, εklm 1⁄4 1 when the sequence fklmg is an even permutation of f123g, and εklm 1⁄4 0 otherwise. The position vector is denoted x 1⁄4 xmim 1⁄4 x1i1 þ x2i2 þ x3i3, where i1; i2; i3 denote the base vectors of the righthanded Cartesian reference frame. Whenever convenient, the vector x is written in subscript notation xk. Time invariance is exploited by using a time-Laplace transformation defined on a space-time dependent vector function Eðx; tÞ as Êðx; sÞ 1⁄4 Z ∞ t1⁄40 Eðx; tÞ expð−stÞdt; (1) where s denotes the Laplace transformation variable, which can be taken real and positive, or s 1⁄4 iω, with ω 1⁄4 2πf, f being the natural frequency. Horizontal shift-invariance, as present in a 1D medium, allows carrying out a 2D spatial Fourier transformation, and it is defined on a Laplace-transformed vector function Êðx; sÞ as ~ EðkT; z; sÞ 1⁄4 Z ∞ xT1⁄4−∞ Êðx; sÞ expð−ikαxαÞdxT; (2) where the subscript T is used to indicate the horizontal vector and kα denotes the two components of the horizontal wavenumber vector. The horizontal and vertical conductivities are denoted as η; ηðvÞ and η 1⁄4 σ þ sε, ηðvÞ 1⁄4 σðvÞ þ sεðvÞ, where σ; σðvÞ are the horizontal and vertical electric conductivity, ε; εðvÞ are the horizontal and vertical electric permittivity. The horizontal and vertical transverse resistivity functions are denoted as ζ; ζðvÞ and ζ 1⁄4 sμ, ζðvÞ 1⁄4 sμðvÞ, with magnetic permeability μ. The wavenumbers related to the vertical conductivity and to the transverse resistivity, are denoted γ; γ̄ with γ 1⁄4 ffiffiffiffiffiffiffiffiffi ζηðvÞ p and γ̄ 1⁄4 ffiffiffiffiffiffiffiffiffi ζðvÞη p . The configuration has an electric or magnetic current type source located in a layer of finite thickness that is embedded in two different layered half-spaces. The medium is a stack of N layers of finite thickness that is bounded at the top and bottom by homogeneous half-spaces. For example, the upper half-space can act as a model for air, and the lower half-space can be taken to lie outside the zone of interest, such that it can act as a radiation boundary condition. This stack of layers is a model for a 1D layered earth. Each layer is VTI characterized by horizontal electric and magnetic parameters that differ from the vertical ones. This type of anisotropy allows complete splitting between TE and TM modes. The solution is obtained in a form that requires the least amount of derivations and the corresponding code requires the least amount of numerical operations. In a horizontally layered VTI earth model the vertical direction is the only direction in which the medium parameters are piecewise constant, while in both horizontal directions they are constant. The vertical direction is therefore chosen as the direction of reference to determine TE modes and TM modes. The TE mode has no vertical electric field component and the TM mode has no vertical magnetic field component. The TE mode and TM mode wave field equations are therefore independently obtained for the vertical components of the magnetic and electric field strengths, respectively. Here, all electromagnetic Green’s functions are derived in the horizontal wavenumber domain for the electric and magnetic fields generated by an electric or magnetic point source somewhere in the layered VTI medium. The receiver is also located somewhere in the VTI medium. Once all Green’s functions are found, the electric field generated by an electric dipole is derived in the space frequency domain by means of an integral equation. Finally, we illustrate the solution method using numerical examples. PROPERTIES OF THE ELECTROMAGNETIC FIELD AND THE REDUCTION OF WORK Maxwell’s equations can be written in the space-Laplace transformed domain as −εkmp∂mĤpðx; sÞ þ ηkrðx; sÞÊrðx; sÞ 1⁄4 −Ĵkðx; sÞ; (3) εjnr∂nÊrðx; sÞ þ ζjpðx; sÞĤpðx; sÞ 1⁄4 −Ĵj ðx; sÞ; (4) which are valid for arbitrary anisotropic media with relaxation and losses. If we substitute Êðx; sÞ 1⁄4 Ĥðx; sÞ, ηmnðx; sÞ 1⁄4 −ζmnðx; sÞ, and Ĵðx; sÞ 1⁄4 −Ĵðx; sÞ, the same Maxwell equations are obtained. This is known as the duality principle. This implies that for any partial or full solution of a boundary value problem that is found for the electric field vector Êðx; sÞ, or a component thereof, the solution for the magnetic field vector Ĥðx; sÞ or the same component thereof can be written down using these substitutions. To reduce bookkeeping, source-receiver reciprocity can be used. Let the Green’s tensor functions related to an electric current point dipole satisfy the following equations: −εkmp∂mĜ me ps ðx; x 0; sÞ þ ηkrðx; sÞĜ rs ðx; x 0; sÞ 1⁄4 −δksðx − x 0Þ; (5) εjnr∂nĜ ee rqðx; x 0; sÞ þ ζjpðx; sÞĜ pq ðx; x 0; sÞ 1⁄4 0; (6) −εkmp∂mĜ mm ps ðx; x 0; sÞ þ ηkrðx; sÞĜ rs ðx; x 0; sÞ 1⁄4 0; (7) F2 Hunziker et al. D ow nl oa de d 04 /2 0/ 15 to 1 31 .1 80 .1 31 .2 42 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / εjnr∂nĜ em rq ðx; x 0; sÞ þ ζjpðx; sÞĜ pq ðx; x 0; sÞ 1⁄4 −δjqðx − x 0Þ; (8) source-receiver reciprocity then states: Ĝ rs ðx; x 0; sÞ 1⁄4 Ĝ sr ðx 0; x; sÞ; (9) Ĝ rq ðx; x 0; sÞ 1⁄4 −Ĝ qr ðx 0; x; sÞ: (10) In these equations the Green’s tensors have superscripts, in which the first superscript denotes the field type and the second superscript denotes the source type. Similarly, the first subscript denotes the field vector component and the second subscript denotes the source vector component. From the duality principle, it was already established that Ĝ pq ðx; x 0; s; ηkr; ζijÞ 1⁄4 −Ĝ pqðx; x 0; s;−ζkr;−ηijÞ; (11) Ĝ pq ðx; x 0; s; ηkr; ζijÞ 1⁄4 −Ĝ pq ðx; x 0; s;−ζkr;−ηijÞ: (12) From these results, it is clear that the solution for the magnetic field Green’s functions are known when those for the electric fields have been found. Moreover, knowing six of the nine components of the Green’s tensor representing the electric field generated by an electric current dipole suffices to know all components of the Green’s tensor describing the electric field generated by an electric source dipole and at the same time all components of the Green’s tensor describing the magnetic field generated by a magnetic source dipole. To break down the problem, first, Maxwell’s equations are transformed to the horizontal wavenumber domain. This simplifies the analysis because Maxwell’s equations are reduced to one scalar second-order differential equation for the vertical electric field component and one for the vertical magnetic field component. This leads to further reduction in work because the whole electromagnetic field is known when the vertical component of the electric field is known for all sources. This is demonstrated in the next section. DECOMPOSITION OF THE ELECTROMAGNETIC FIELD AND THE FREQUENCY WAVENUMBER DOMAIN HOMOGENEOUS SPACE GREEN’S FUNCTIONS We start with the solution in a homogeneous space, but we use the coupled Maxwell equations in the horizontal wavenumber domain to prepare for the layered medium treated afterward. They are given by ηðvÞ ~ E3 1⁄4 − ~ J3 − ε3λβikλ ~ Hβ; (13) ζðvÞ ~ H3 1⁄4 − ~ J3 þ ε3λβikλ ~ Eβ; (14) η ~ Eα 1⁄4 − ~ Jα þ εα3β∂3 ~ Hβ − εαλ3ikλ ~ H3; (15) ζ ~ Hβ 1⁄4 − ~ Jβ − εβ3α∂3 ~ Eα þ εβλ3ikλ ~ E3: (16) The horizontal electric field components are eliminated by substituting the expression of equation 15 into 16: 1⁄2ηζ − ∂3∂3 ~ Hα 1⁄4 −η ~ Jα þ εα3β∂3 ~ Jβ þ ikα∂3 ~ H3 þ ηεαλ3ikλ ~ E3; (17) and using equation 17 in equation 13, to obtain ðηζ − ∂3∂3ÞηðvÞ ~ E3 1⁄4 ð∂3∂3 − ηζÞ ~ J3 þ ηε3αβikα ~ Jβ − ε3αβikαðεβ3α∂3 ~ Jα þ ikβ∂3 ~ H3 þ ηεβλ3ikλ ~ E3Þ; (18) leading to the following equation for the vertical electric field component ð∂3∂3 − Γ2ÞηðvÞ ~ E3 1⁄4 ηζ ~ J3 þ ∂3ðikα ~ Jα − ∂3 ~ J3Þ − ηε3αβikα ~ Jβ ; (19) with vertical wavenumber related to vertical conductivity given by Γ 1⁄4 ffiffiffiffiffiffiffiffiffiffiffiffi η∕ηðvÞ p ffiffiffiffiffiffiffiffiffiffiffiffiffiffi κ þ γ p , with κ 1⁄4 ðk1 þ k2Þ being the horizontal radial wavenumber, and RfΓg ≥ 0. To arrive at equations 17 and 19, the epsilon-delta identity has been used, εkmpεpnr 1⁄4 δknδmr − δkrδmn: (20) Without loss of generality, we take the sources to be point sources, defined by f ~ Jk; ~ Jk g 1⁄4 fÎkðsÞ; Îk ðsÞgδðx3 − x3Þ; (21) with the source position vector given by x3 1⁄4 f0; 0; zSg and the frequency spectrum of the source specified by Î k ðsÞ. The two corresponding Green’s functions are solutions of modified Helmholtz equations, ð∂3∂3 − Γ2Þ ~ G 1⁄4 −δðx3 − x3Þ: (22) The solutions to these equations are well known, and given by ~ Gðx3 − x3Þ 1⁄4 expð−ΓhÞ 2Γ ; (23) where h 1⁄4 jx3 − x3 j. Equation 23 is the TM-mode Green’s function. In view of equations 19, 21, and 23, the vertical component of the electric field is given by ηðvÞ ~ E3 1⁄4 ð−ηζĴ3 − ∂31⁄2ikαĴα − ∂3Ĵ3 þ ηε3αβikαĴβ Þ ~ Gðx3 − x3Þ: (24) Working out the derivatives gives ~ E3 1⁄4 X 3 ~ Gðx3 − x3Þ − ~ J3 ηðvÞ ; (25) with the source factor given by The EM response in a 1D VTI medium F3 D ow nl oa de d 04 /2 0/ 15 to 1 31 .1 80 .1 31 .2 42 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / X 3 1⁄4 Γ ηðvÞ signðx3 − x3ÞikαĴα þ ηκ ðηðvÞÞ2 Ĵ e 3 þ η ηðvÞ ε3αβikαĴ m β ; (26) where the sign function is defined as signðx3 − x3 0Þ 1⁄4 8< : −1 for x3 < x3; 0 for x3 1⁄4 x3; 1 for x3 > x3: The vertical component of the magnetic field is found by direct substitution as described above. To obtain the horizontal electric field components from the vertical electric and magnetic field components, first, both sides of equation 15 are divided by η and multiplied by ikλikα ikλikα ~ Eα 1⁄4 −ηikλ1⁄2ikα ~ Jα − ikαεα3β∂3 ~ Hβ ; (27) and then equation 13 is used to eliminate the horizontal components of the magnetic field ikλikα ~ Eα 1⁄4 −ηikλ1⁄2ikα ~ Jα − ∂3ðηðvÞ ~ E3 þ ~ J3Þ : (28) Then, both sides of equation 14 are multiplied by εαβ3ikβ to obtain ikαikλ ~ Eλ − ikλikλ ~ Eα 1⁄4 εαβ3ikβðζðvÞ ~ H3 þ ~ J3 Þ: (29) Next, equation 28 is used to eliminate ikαikλ ~ Eλ from equation 29: ~ Eα 1⁄4 ikα ηκ 1⁄2ikβ ~ Jβ − ∂3ðηðvÞ ~ E3 þ ~ J3Þ þ εαβ3 ikβ κ ðζðvÞ ~ H3 þ ~ J3 Þ: (30) In equation 30, the TM mode and the TE mode are separated from each other because ~ E3 is a pure TM mode and ~ H3 is a pure TE mode. With the aid of equations 25 and 26, and the corresponding equations for the vertical component of the magnetic field that are obtained from the duality principle, the horizontal electric field components of equation 30 can be written as ~ Eα 1⁄4 X α ~ Gþ X α ~̄ G; (31) where the TE mode Green’s function is given by ~̄ Gðx3 − x3Þ 1⁄4 expð−Γ̄hÞ 2Γ̄ ; (32) with the vertical wavenumber given by Γ̄ 1⁄4 ffiffiffiffiffiffiffiffiffiffiffiffi ζ∕ζðvÞ p ffiffiffiffiffiffiffiffiffiffiffiffiffiffi κ þ γ̄ p . The electric-field-related TM mode and TE mode horizontal source components are given by X α 1⁄4 ikαikβΓĴ e β ηκ þ ikαΓĴ e 3 signðx3 − x3Þ ηðvÞ þ ikαε3λβikλΓĴ m β signðx3 − x3Þ κ ; (33) X α 1⁄4 −ζ ikαikβĴ e β þ κĴα κ þ εαλ3ikλikβΓ̄Ĵ m β signðx3 − x3Þ κ þ ζ ζðvÞ εαβ3ikβĴ m 3 : (34) With the above results the Green’s tensors are found as ~ G kr 1⁄4 0 BBB@ ðik1ÞΓ ηκ ik1ik2Γ ηκ − ik1∂3 ηðvÞ ik1ik2Γ ηκ ðik2ÞΓ ηκ − ik2∂3 ηðvÞ − ik1∂3 ηðvÞ − ik2∂3 ηðvÞ ∂3∂3−ηζ ηðvÞ 1 CCCA ~ G þ 0 BB@ ζðik2Þ κ − ζik1ik2 κ 0 − ζik1ik2 κ ζðik1Þ κ 0 0 0 0 1 CCA ~̄ G; (35) ~ G kr 1⁄4 0 BB@ ik1ik2∂3 κ − ðik1Þ ∂3 κ 0 ðik2Þ∂3 κ − ik1ik2∂3 κ 0 − ηik2 ηðvÞ ηik1 ηðvÞ 0 1 CCA ~ G þ 0 BB@ − ik1ik2∂3 κ − ðik2Þ ∂3 κ ζik2 ζðvÞ ðik1Þ∂3 κ ik1ik2∂3 κ − ζik1 ζðvÞ 0 0 0 1 CCA ~̄ G: (36) This decomposition demonstrates that once the vertical electric and magnetic field components are known, the entire electromagnetic field is known. From the duality principle, it is clear that if the vertical electric field is known, the vertical magnetic field is also known. The same applies to the horizontal components of the magnetic fields. Hence, finding the solution for the vertical electric field in the layered earth model suffices to know the entire electromagnetic field. THE SPACE FREQUENCY DOMAIN HOMOGENEOUS SPACE GREEN’S FUNCTIONS The elements of equations 35 and 36 correspond to homogeneous space Green’s functions and can be transformed back to space-frequency in closed form. This is done in this section and results in similar expressions as, for example, found by Weiglhofer (1990) or Abubakar and Habashy (2006). In the first element of ~ G αβ, the factor Γ∕ðηκ2Þ is rewritten as Γ ηκ 1⁄4 1 ηðvÞΓ þ ζ κΓ : (37) With this substitution the Green’s function corresponding to the electric field generated by an electric current source in space domain can be written as an inverse spatial Fourier-Bessel transformation F4 Hunziker et al. D ow nl oa de d 04 /2 0/ 15 to 1 31 .1 80 .1 31 .2 42 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / Ĝ αβðx; xS;ωÞ 1⁄4 ∂α∂β 4πηðvÞ Z ∞ κ1⁄40 expð−ΓhÞ Γ J0ðκrÞκdκ þ ζ∂α∂β 4π Z ∞ κ1⁄40 expð−ΓhÞ Γ − expð−Γ̄hÞ Γ̄ J0ðκrÞκdκ − ζδαβ 4π Z ∞ κ1⁄40 expð−Γ̄hÞ Γ̄ J0ðκrÞκdκ; (38) where r 1⁄4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x þ y p is the horizontal distance between source and receiver and J0 denotes the zero-order Bessel function. The first and third integrals are standard Fourier-Bessel transforms, and the second can be evaluated by first performing one derivative on the Bessel function in the integrands. That leads to Ĝðx; xSÞ 1⁄4 1 4π Z ∞ κ1⁄40 expð−ΓhÞ Γ J0ðκrÞκdκ 1⁄4 expð−γRÞ 4πR ffiffiffiffiffiffiffiffiffiffiffiffi η∕ηðvÞ p ;